Questionnaire on Well-Being (QWB)

Reproducing EFA

Author
Affiliation

Magnus Johansson

RISE Research Institutes of Sweden

Published

2024-10-04

Code
library(readxl)
library(tidyverse)
library(RISEkbmRasch) # devtools::install_github("pgmj/RISEkbmRasch")

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename

Questionnaire on Well-Being (QWB), 18 items, each item is scored on a scale of 0 to 4 [@hlynsson_evaluating_2024]. Data from the same paper. We’ll use the data from the first study that was used in the EFA in the paper.

Code
### import data - this is just sample code, the files do not exist
d1 <- read_excel("data/study_one.xlsx") # replace with your datafile as needed

### Load item information
# make sure that variable names in df match with itemlabels$itemnr
itemlabels <- data.frame(
  stringsAsFactors = FALSE,
  itemnr = paste0("swb", 1:18),
  item = c(
    "... felt calm and relaxed?",
    "... had a good appetite and enjoyed eating?",
    "... have been able to take initiatives and get started with what you wanted to do?",
    "... felt optimistic and viewed things on the bright side?",
    "... felt active and filled with energy?",
    "... felt a strong zest for life?",
    "... been able to object and to assert yourself when this is needed?",
    "... felt satisfied with your life in its present situation?",
    "... felt that your life is meaningful?",
    "... slept well and got the right amount of sleep?",
    "... have been able to be focused and concentrated on today’s tasks?",
    "... felt interested in various activities and in people around you?",
    "... felt happy and harmonious?",
    "... felt satisfied with yourself?",
    "... been able to make decisions and carry them out?",
    "... been able to stay in the here and now and to let go of thoughts about problems?",
    "... had the power to recover one’s strength if something has been stressful or difficult?",
    "... felt that you are well and healthy?"
  )
)
#write_csv(itemlabels,"data/itemlabels_swb.csv")
Code
# swb items extract
df <- d1 %>% 
  select(starts_with("qwb")) %>% 
  select(contains("screening")) %>% 
  select(!contains("sum")) %>% 
  set_names(itemlabels$itemnr)
Code

Code
RImissingP(df,n = 1000)

All participants that have missing data are completely missing (no missing data on item level), except one respondent with 33.3% missing. We will remove all those with any missing responses from our data.

Code
nrow(df)
[1] 659
Code 
[1] 494
Code
df <- na.omit(df)

We end up with 494 respondents, compared to the 495 in the paper, which seems likely due to them keeping the one person with missing item responses in the dataset. This should not contribute to any noticable differences in analysis outcomes.

1 All items in the analysis

Code
itemnr item
swb1 ... felt calm and relaxed?
swb2 ... had a good appetite and enjoyed eating?
swb3 ... have been able to take initiatives and get started with what you wanted to do?
swb4 ... felt optimistic and viewed things on the bright side?
swb5 ... felt active and filled with energy?
swb6 ... felt a strong zest for life?
swb7 ... been able to object and to assert yourself when this is needed?
swb8 ... felt satisfied with your life in its present situation?
swb9 ... felt that your life is meaningful?
swb10 ... slept well and got the right amount of sleep?
swb11 ... have been able to be focused and concentrated on today’s tasks?
swb12 ... felt interested in various activities and in people around you?
swb13 ... felt happy and harmonious?
swb14 ... felt satisfied with yourself?
swb15 ... been able to make decisions and carry them out?
swb16 ... been able to stay in the here and now and to let go of thoughts about problems?
swb17 ... had the power to recover one’s strength if something has been stressful or difficult?
swb18 ... felt that you are well and healthy?

2 Descriptives of raw data

Response distribution for all items are summarized below.

Code
Response category Number of responses Percent
0 537 6.0
1 2432 27.4
2 3422 38.5
3 1981 22.3
4 520 5.8

2.1 Descriptives - item level

Code
RIlistItemsMargin(df, fontsize = 12)
itemnr item
swb1 … felt calm and relaxed?
swb2 … had a good appetite and enjoyed eating?
swb3 … have been able to take initiatives and get started with what you wanted to do?
swb4 … felt optimistic and viewed things on the bright side?
swb5 … felt active and filled with energy?
swb6 … felt a strong zest for life?
swb7 … been able to object and to assert yourself when this is needed?
swb8 … felt satisfied with your life in its present situation?
swb9 … felt that your life is meaningful?
swb10 … slept well and got the right amount of sleep?
swb11 … have been able to be focused and concentrated on today’s tasks?
swb12 … felt interested in various activities and in people around you?
swb13 … felt happy and harmonious?
swb14 … felt satisfied with yourself?
swb15 … been able to make decisions and carry them out?
swb16 … been able to stay in the here and now and to let go of thoughts about problems?
swb17 … had the power to recover one’s strength if something has been stressful or difficult?
swb18 … felt that you are well and healthy?
Code

Code

Code

Code
library(TAM)
tam1 <- tam(as.matrix(df), irtmodel = "PCM", verbose = FALSE) # run TAM Rasch Partial Credit Model on our data, which uses Marginal Maximum Likelihood estimation
plot(tam1) # create ICC plots
Iteration in WLE/MLE estimation  1   | Maximal change  1.5065 
Iteration in WLE/MLE estimation  2   | Maximal change  0.8894 
Iteration in WLE/MLE estimation  3   | Maximal change  0.4824 
Iteration in WLE/MLE estimation  4   | Maximal change  0.1072 
Iteration in WLE/MLE estimation  5   | Maximal change  0.0011 
Iteration in WLE/MLE estimation  6   | Maximal change  0.0001 
----
 WLE Reliability= 0.933 

....................................................
 Plots exported in png format into folder:
 /Users/magnuspjo/Documents/Hub/rasch_välmåendeskalan/Plots

Some cells have few (4-6) responses. No items need reverse scoring.

3 Exploratory factor analysis

Let’s look at the data using EFA. A lot of the code for this analysis was borrowed from https://solomonkurz.netlify.app/blog/2021-05-11-yes-you-can-fit-an-exploratory-factor-analysis-with-lavaan/

Since the paper used varimax rotation and WLSMV estimator, we will do the same.

Code
library(lavaan)
f1 <- 'efa("efa")*f1 =~ swb1 + swb2 + swb3 + swb4 + swb5 + swb6 + swb7 + swb8 +
                    swb9 + swb10 + swb11 + swb12 + swb13 + swb14 + swb15 + 
                    swb16 + swb17 + swb18'

# 2-factor model
f2 <- 'efa("efa")*f1 + 
       efa("efa")*f2 =~ swb1 + swb2 + swb3 + swb4 + swb5 + swb6 + swb7 + swb8 +
                    swb9 + swb10 + swb11 + swb12 + swb13 + swb14 + swb15 + 
                    swb16 + swb17 + swb18'

# 3-factor
f3 <- '
efa("efa")*f1 +
efa("efa")*f2 +
efa("efa")*f3 =~ swb1 + swb2 + swb3 + swb4 + swb5 + swb6 + swb7 + swb8 +
                    swb9 + swb10 + swb11 + swb12 + swb13 + swb14 + swb15 + 
                    swb16 + swb17 + swb18'

# 4-factor
f4 <- '
efa("efa")*f1 +
efa("efa")*f2 +
efa("efa")*f3 +
efa("efa")*f4 =~ swb1 + swb2 + swb3 + swb4 + swb5 + swb6 + swb7 + swb8 +
                    swb9 + swb10 + swb11 + swb12 + swb13 + swb14 + swb15 + 
                    swb16 + swb17 + swb18'

# 5-factor
f5 <- '
efa("efa")*f1 +
efa("efa")*f2 +
efa("efa")*f3 +
efa("efa")*f4 + 
efa("efa")*f5 =~ swb1 + swb2 + swb3 + swb4 + swb5 + swb6 + swb7 + swb8 +
                    swb9 + swb10 + swb11 + swb12 + swb13 + swb14 + swb15 + 
                    swb16 + swb17 + swb18'

efa_f1 <- 
  cfa(model = f1,
      data = df,
      rotation = "varimax",
      estimator = "WLSMV",
      ordered = TRUE)
efa_f2 <- 
  cfa(model = f2,
      data = df,
      rotation = "varimax",
      estimator = "WLSMV",
      ordered = TRUE)
efa_f3 <- 
  cfa(model = f3,
      data = df,
      rotation = "varimax",
      estimator = "WLSMV",
      ordered = TRUE)
efa_f4 <- 
  cfa(model = f4,
      data = df,
      rotation = "varimax",
      estimator = "WLSMV",
      ordered = TRUE)
efa_f5 <- 
  cfa(model = f5,
      data = df,
      rotation = "varimax",
      estimator = "WLSMV",
      ordered = TRUE)

3.1 Scaled fit metrics

For WLSMV, the .scaled metrics should be reported.

Code
fit_metrics_scaled <- c("chisq.scaled", "df", "pvalue.scaled", 
                        "cfi.scaled", "tli.scaled", "rmsea.scaled", 
                        "rmsea.ci.lower.scaled","rmsea.ci.upper.scaled",
                        "srmr")


rbind(
  fitmeasures(efa_f1, fit_metrics_scaled),
  fitmeasures(efa_f2, fit_metrics_scaled),
  fitmeasures(efa_f3, fit_metrics_scaled),
  fitmeasures(efa_f4, fit_metrics_scaled),
  fitmeasures(efa_f5, fit_metrics_scaled)
  ) %>% 
  as.data.frame() %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  rename(Chi2.scaled = chisq.scaled,
         p.scaled = pvalue.scaled,
         CFI.scaled = cfi.scaled,
         TLI.scaled = tli.scaled,
         RMSEA.scaled = rmsea.scaled,
         CI_low.scaled = rmsea.ci.lower.scaled,
         CI_high.scaled = rmsea.ci.upper.scaled,
         SRMR = srmr) %>% 
  add_column(Model = paste0(1:5,"-factor"), .before = "Chi2.scaled") %>% 
  knitr::kable()
Model Chi2.scaled df p.scaled CFI.scaled TLI.scaled RMSEA.scaled CI_low.scaled CI_high.scaled SRMR
1-factor 894.634 135 0 0.944 0.937 0.107 0.100 0.114 0.061
2-factor 642.016 118 0 0.961 0.950 0.095 0.088 0.102 0.050
3-factor 390.419 102 0 0.979 0.968 0.076 0.068 0.084 0.035
4-factor 255.405 87 0 0.988 0.978 0.063 0.054 0.072 0.028
5-factor 172.833 73 0 0.993 0.985 0.053 0.043 0.063 0.021

Since these are ordinal data used with WLSMV, we cannot apply rule of thumb cutoff values such as Hu & Bentler [-@hu_cutoff_1999], since those are based on continuous data and maximum likelihood estimation. But we can see what happens with the different metrics in the different models and how they compare.

3.2 Factor loadings 1-factor

Code
standardizedsolution(efa_f1) %>% 
  filter(op == "=~") %>% 
  mutate(item  = str_remove(rhs, "swb") %>% as.double(),
         factor = str_remove(lhs, "f"))
   lhs op   rhs est.std    se      z pvalue ci.lower ci.upper item factor
1   f1 =~  swb1   0.783 0.019 40.224      0    0.745    0.821    1      1
2   f1 =~  swb2   0.566 0.031 18.301      0    0.505    0.627    2      1
3   f1 =~  swb3   0.534 0.031 17.199      0    0.473    0.595    3      1
4   f1 =~  swb4   0.715 0.024 30.156      0    0.668    0.761    4      1
5   f1 =~  swb5   0.740 0.022 34.159      0    0.698    0.783    5      1
6   f1 =~  swb6   0.743 0.022 34.237      0    0.700    0.785    6      1
7   f1 =~  swb7   0.815 0.016 52.008      0    0.784    0.846    7      1
8   f1 =~  swb8   0.873 0.014 63.481      0    0.846    0.900    8      1
9   f1 =~  swb9   0.784 0.021 37.401      0    0.743    0.825    9      1
10  f1 =~ swb10   0.784 0.020 39.219      0    0.745    0.823   10      1
11  f1 =~ swb11   0.778 0.020 38.898      0    0.739    0.817   11      1
12  f1 =~ swb12   0.679 0.026 25.989      0    0.628    0.730   12      1
13  f1 =~ swb13   0.404 0.039 10.468      0    0.328    0.479   13      1
14  f1 =~ swb14   0.627 0.029 21.782      0    0.571    0.683   14      1
15  f1 =~ swb15   0.777 0.018 42.393      0    0.742    0.813   15      1
16  f1 =~ swb16   0.687 0.026 26.849      0    0.637    0.737   16      1
17  f1 =~ swb17   0.769 0.019 40.441      0    0.732    0.807   17      1
18  f1 =~ swb18   0.685 0.024 28.688      0    0.638    0.731   18      1

All loadings (est.std in the table above) are above 0.4 in the one-factor model. Let’s review the modification indices, filtered to only include those with mi/chi2 values above 15.

3.3 Modification indices 1-factor

Code
modificationIndices(efa_f1,
                    standardized = T) %>% 
  as.data.frame(row.names = NULL) %>% 
  filter(mi > 15) %>% 
  arrange(desc(mi)) %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  knitr::kable()
lhs op rhs mi epc sepc.lv sepc.all sepc.nox
swb4 ~~ swb5 64.432 -0.210 -0.210 -0.448 -0.448
swb7 ~~ swb8 47.569 -0.161 -0.161 -0.568 -0.568
swb1 ~~ swb2 34.092 -0.190 -0.190 -0.371 -0.371
swb15 ~~ swb17 27.174 -0.140 -0.140 -0.348 -0.348
swb5 ~~ swb12 22.817 -0.138 -0.138 -0.281 -0.281
swb11 ~~ swb17 22.577 -0.131 -0.131 -0.326 -0.326
swb1 ~~ swb16 20.692 -0.142 -0.142 -0.313 -0.313
swb5 ~~ swb6 19.294 -0.123 -0.123 -0.273 -0.273
swb1 ~~ swb14 17.851 -0.155 -0.155 -0.320 -0.320
swb7 ~~ swb11 15.753 -0.102 -0.102 -0.280 -0.280
swb10 ~~ swb15 15.145 -0.109 -0.109 -0.279 -0.279

Here we can see several large residual correlations. These are the top 8:

  • items 4 and 5
  • 7 and 8
  • 1 and 2
  • 15 and 17
  • 5 and 12
  • 11 and 17
  • 1 and 16
  • 5 and 6

There seems to be some clustering going on here, which indicates multiple dimensions.

Let’s look at the 3 factor model as a comparison, and review how the items listed above fit in the 3 factors.

3.4 Plot 3-factor model

In the paper, the lowest standardized factor loading considered acceptable is 0.3, so we will use that here as well.

Code
standardizedsolution(efa_f3) %>% 
  filter(op == "=~") %>% 
  mutate(item  = str_remove(rhs, "swb") %>% as.double(),
         factor = str_remove(lhs, "f")) %>% 
  # plot
  ggplot(aes(x = est.std, xmin = ci.lower, xmax = ci.upper, y = item)) +
  annotate(geom = "rect",
           xmin = -1, xmax = 1,
           ymin = -Inf, ymax = Inf,
           fill = "grey90") +
  annotate(geom = "rect",
           xmin = -0.7, xmax = 0.7,
           ymin = -Inf, ymax = Inf,
           fill = "grey93") +
  annotate(geom = "rect",
           xmin = -0.3, xmax = 0.3,
           ymin = -Inf, ymax = Inf,
           fill = "grey96") +
  geom_vline(xintercept = 0, color = "white") +
  geom_pointrange(aes(alpha = abs(est.std) < 0.3),
                  fatten = 10) +
  geom_text(aes(label = item, color = abs(est.std) < 0.3),
            size = 4) +
  scale_color_manual(values = c("white", "transparent")) +
  scale_alpha_manual(values = c(1, 1/3)) +
  scale_x_continuous(expression(lambda[standardized]), 
                     expand = c(0, 0), limits = c(-1, 1),
                     breaks = c(-1, -0.7, -0.3, 0, 0.3, 0.7, 1),
                     labels = c("-1", "-.7", "-.3", "0", ".3", ".7", "1")) +
  scale_y_continuous(breaks = 1:18, sec.axis = sec_axis(~ . * 1, breaks = 1:18)) +
  ggtitle("Factor loadings for the 3-factor model") +
  theme(legend.position = "none") +
  facet_wrap(~ factor, labeller = label_both)

Several crossloadings are apparent when using the low cutoff of 0.3. Let’s raise the bar to 0.4 and perhaps achieve more clarity.

Code
standardizedsolution(efa_f3) %>% 
  filter(op == "=~") %>% 
  mutate(item  = str_remove(rhs, "swb") %>% as.double(),
         factor = str_remove(lhs, "f")) %>% 
  # plot
  ggplot(aes(x = est.std, xmin = ci.lower, xmax = ci.upper, y = item)) +
  annotate(geom = "rect",
           xmin = -1, xmax = 1,
           ymin = -Inf, ymax = Inf,
           fill = "grey90") +
  annotate(geom = "rect",
           xmin = -0.7, xmax = 0.7,
           ymin = -Inf, ymax = Inf,
           fill = "grey93") +
  annotate(geom = "rect",
           xmin = -0.4, xmax = 0.4,
           ymin = -Inf, ymax = Inf,
           fill = "grey96") +
  geom_vline(xintercept = 0, color = "white") +
  geom_pointrange(aes(alpha = abs(est.std) < 0.4),
                  fatten = 10) +
  geom_text(aes(label = item, color = abs(est.std) < 0.4),
            size = 4) +
  scale_color_manual(values = c("white", "transparent")) +
  scale_alpha_manual(values = c(1, 1/3)) +
  scale_x_continuous(expression(lambda[standardized]), 
                     expand = c(0, 0), limits = c(-1, 1),
                     breaks = c(-1, -0.7, -0.4, 0, 0.4, 0.7, 1),
                     labels = c("-1", "-.7", "-.4", "0", ".4", ".7", "1")) +
  scale_y_continuous(breaks = 1:18, sec.axis = sec_axis(~ . * 1, breaks = 1:18)) +
  ggtitle("Factor loadings for the 3-factor model") +
  theme(legend.position = "none") +
  facet_wrap(~ factor, labeller = label_both)

These are the residual correlations from the 1-factor model sorted into factors in the model above:

factor 1:

  • 4 and 5
  • 5 and 6
  • 5 and 12

factor 2:

  • 1 and 2
  • 1 and 16

factor 3:

  • 7 and 8
  • 15 and 17
  • 11 and 17

Item pairs with residual correlations indicates that either the items make up a separate dimension/cluster in data, and/or they are too similar in content and only one of the two items can be included in the scale.

For unidimensionality to hold, there cannot be local dependence between items - items should only be related through the latent variable. Residuals are the variance in an item that is not explained by the latent variable, and thus the pattern of residuals is important to assess local dependence. Residual correlations are such a pattern.

3.5 Modification indices 3-factor

Code
modificationIndices(efa_f3,
                    standardized = T) %>% 
  as.data.frame(row.names = NULL) %>% 
  filter(mi > 5) %>% 
  arrange(desc(mi)) %>% 
  mutate(across(where(is.numeric),~ round(.x, 3))) %>%
  knitr::kable()
lhs op rhs mi epc sepc.lv sepc.all sepc.nox
swb15 ~~ swb17 24.233 -0.144 -0.144 -0.389 -0.389
swb7 ~~ swb8 16.739 -0.148 -0.148 -0.729 -0.729
swb6 ~~ swb15 14.164 0.145 0.145 0.371 0.371
swb7 ~~ swb15 12.128 0.118 0.118 0.380 0.380
swb6 ~~ swb7 10.764 -0.114 -0.114 -0.355 -0.355
swb10 ~~ swb15 10.318 -0.095 -0.095 -0.258 -0.258
swb12 ~~ swb13 9.720 -0.146 -0.146 -0.233 -0.233
swb8 ~~ swb17 7.809 0.104 0.104 0.430 0.430
swb5 ~~ swb6 7.246 -0.113 -0.113 -0.347 -0.347
swb9 ~~ swb14 7.198 0.114 0.114 0.250 0.250
swb7 ~~ swb17 7.063 0.095 0.095 0.313 0.313
swb5 ~~ swb18 6.130 0.117 0.117 0.327 0.327
swb2 ~~ swb9 5.476 -0.093 -0.093 -0.208 -0.208
swb16 ~~ swb18 5.095 -0.088 -0.088 -0.191 -0.191

A lot of residual correlations, particularly in the 3rd factor, which also has the most items.

4 Summary comments

These 18 items are clearly not unidimensional and should not be sum scored. Further work is needed to assess dimensionality to arrive at a set of items that fulfill unidimensionality requirements.

5 Software used

Code 
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Stockholm
tzcode source: internal

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] lavaan_0.6-18        TAM_4.2-21           CDM_8.2-6           
 [4] RISEkbmRasch_0.2.4.4 janitor_2.2.0        iarm_0.4.3          
 [7] hexbin_1.28.4        catR_3.17            glue_1.7.0          
[10] ggrepel_0.9.6        patchwork_1.3.0      reshape_0.8.9       
[13] matrixStats_1.4.1    psychotree_0.16-1    psychotools_0.7-4   
[16] partykit_1.2-22      mvtnorm_1.3-1        libcoin_1.0-10      
[19] psych_2.4.6.26       mirt_1.42            lattice_0.22-6      
[22] eRm_1.0-6            kableExtra_1.4.0     formattable_0.2.1   
[25] lubridate_1.9.3      forcats_1.0.0        stringr_1.5.1       
[28] dplyr_1.1.4          purrr_1.0.2          readr_2.1.5         
[31] tidyr_1.3.1          tibble_3.2.1         ggplot2_3.5.1       
[34] tidyverse_2.0.0      readxl_1.4.3        

loaded via a namespace (and not attached):
  [1] vcd_1.4-12           rstudioapi_0.16.0    audio_0.1-11        
  [4] jsonlite_1.8.9       magrittr_2.0.3       farver_2.1.2        
  [7] rmarkdown_2.28       vctrs_0.6.5          base64enc_0.1-3     
 [10] htmltools_0.5.8.1    curl_5.2.3           cellranger_1.1.0    
 [13] Formula_1.2-5        dcurver_0.9.2        parallelly_1.38.0   
 [16] htmlwidgets_1.6.4    plyr_1.8.9           testthat_3.2.1.1    
 [19] zoo_1.8-12           admisc_0.35          lifecycle_1.0.4     
 [22] pkgconfig_2.0.3      Matrix_1.7-0         R6_2.5.1            
 [25] fastmap_1.2.0        snakecase_0.11.1     future_1.34.0       
 [28] numDeriv_2016.8-1.1  digest_0.6.37        colorspace_2.1-1    
 [31] rprojroot_2.0.4      Hmisc_5.1-3          vegan_2.6-8         
 [34] labeling_0.4.3       progressr_0.14.0     fansi_1.0.6         
 [37] timechange_0.3.0     abind_1.4-5          mgcv_1.9-1          
 [40] compiler_4.4.1       here_1.0.1           vcdExtra_0.8-5      
 [43] withr_3.0.1          htmlTable_2.4.3      backports_1.5.0     
 [46] carData_3.0-5        relimp_1.0-5         R.utils_2.12.3      
 [49] MASS_7.3-61          sessioninfo_1.2.2    GPArotation_2024.3-1
 [52] permute_0.9-7        pbivnorm_0.6.0       tools_4.4.1         
 [55] foreign_0.8-87       lmtest_0.9-40        future.apply_1.11.2 
 [58] nnet_7.3-19          quadprog_1.5-8       R.oo_1.26.0         
 [61] nlme_3.1-165         inum_1.0-5           checkmate_2.3.2     
 [64] cluster_2.1.6        generics_0.1.3       snow_0.4-4          
 [67] gtable_0.3.5         RPushbullet_0.3.4    tzdb_0.4.0          
 [70] ca_0.71.1            R.methodsS3_1.8.2    data.table_1.16.0   
 [73] hms_1.1.3            car_3.1-2            xml2_1.3.6          
 [76] Deriv_4.1.3          utf8_1.2.4           pillar_1.9.0        
 [79] splines_4.4.1        survival_3.7-0       tidyselect_1.2.1    
 [82] pbapply_1.7-2        knitr_1.48           gridExtra_2.3       
 [85] svglite_2.1.3        xfun_0.46            gnm_1.1-5           
 [88] brio_1.1.5           stringi_1.8.4        yaml_2.3.10         
 [91] evaluate_0.24.0      codetools_0.2-20     beepr_2.0           
 [94] cli_3.6.3            rpart_4.1.23         qvcalc_1.0.3        
 [97] systemfonts_1.1.0    munsell_0.5.1        Rcpp_1.0.13         
[100] globals_0.16.3       polycor_0.8-1        parallel_4.4.1      
[103] listenv_0.9.1        viridisLite_0.4.2    SimDesign_2.17.1    
[106] scales_1.3.0         rlang_1.1.4          mnormt_2.1.1

6 References